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SYSTEM AND METHOD FOR ENABLING SIMULTANEOUS CALIBRATION 

AND IMAGING OF A MEDIUM 

CROSS-REFERENCE TO A RELATED APPLICATION 
[0001] This application claims priority from the Provisional Appln. No. 

60/261,492 filed on January 12, 2001, the entire disclosure of which is incorporated 
herein by reference. 

STATEMENT AS TO FEDERALLY SPONSORED RESEARCH 
[0002] This invention was made with Government support under grants fi-om the 

National Institutes of Health, NIH R29-NS38842 and NIH P41-RR14075, and the U.S. 
Army, No. DAMD 17-99-2-9001. As a result, the Government of the United States of 
America has certain rights to the invention. 

TECHNICAL FIELD 
[0003] The invention relates to imaging, and more particularly to calibration of 

optical methods of imaging highly scattering media. 

BACKGROUND 

[0004] Opaque media such as paint, milk, foam, emulsions, colloidal suspensions 

and human tissue do not strongly absorb visible light, while being optically turbid. This 
turbidity is a result of a very short scattering length characteristic of photons traveling 
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within these media. The light does not travel in straight lines through such substances, 
and instead "diffuses" in a manner similar to heat flow. In other words, photons are 
multiply scattered within these media until they are either absorbed, or they exit the 
boundaries of the medium. 

[0005] Recently, there has been a significant interest in the use of optical 

radiation for imaging within the highly scattering media, such as a biological tissue. 
Photons can travel within a highly scattering medium along a distribution of paths, of 
which very few are straight. Thus, the direction of the light into a highly scattering 
medium and a subsequent direction of the diffuse Ught emitted fi"om the medium provides 
certain information regarding local variations in the scattering and absorption 
coefficients. Such information can identify, for example, a breast or brain tumor, 
bleeding in the brain, or aneurysm. Furthermore, multiple wavelengths can be used to 
spectroscopically determine local tissue concentrations of oxy-hemoglobin (HbO) and 
deoxy-hemoglobin (Hb) in tissue, which can vary in response to some stimulus, e.g., a 
drug. For a general description of such appUcations, see, e.g., A. Yodh et al.. Physics 
Today, 34-40, March 1995. 

[0006] If the spatially varying optical properties of a highly scattering medium 

are known, photon propagation within the medium can be calculated numerically. This 
numerical calculation is simplified when the scattering is much larger than the 
absorption, in which case photon propagation can be approximated as a diffusion process. 
This condition is typically satisfied in a biological tissue in the spectral range of about 

-2- 

NY02:364961.1 




A34927- 069225.0110 



700 to 900 nm. The numerical calculation gives the distribution of hght inside a highly 
scattering medium, and is usually referred to as the "forward calculation." For a medium 
being imaged, however, the "inverse calculation" should be solved, e.g., by deducing the 
distribution of optical properties within the medium from the diffuse light measurements. 
Numerical techniques for performing the inverse calculation include singular value 
decomposition (SVD), simultaneous iterative reconstruction technique (SIRT), K-space 
diffraction tomography, and a use of an extended Kaman filter. For a general review of 
techniques for the forward and inverse calculations, see, e.g., S.R. Arridge, Inverse 
Problems, 15:S41, 1999. 

[0007] In the art of diffuse optical tomography ("DOT"), multiple sources 

sequentially direct the light into a highly scattering medium (e.g., tissue), at spatially 
separated locations. For each such source, multiple detectors on the tissue measure the 
diffuse hght emitted from the sample at spatially separated locations. The detectors may 
further obtain one or more parameters of the diffuse light emitted, e.g., fluence, and then 
utilize those parameters as input in the inverse calculation. However, the measurements 
can include various errors caused by, for example, source and detector coupling to the 
tissue, soxirce and detector positional uncertainties, fluctuations in the source power, and 
variations in the detector gain. 

[0008] to minimize these uncertainties, DOT systems are typically calibrated with 

initial measurements for a known sample, and the calibration is used to correct the 
subsequent measurements for imaging an unknown sample. Unfortunately, errors can 
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vary from a measurement to another measurement because of, e.g., the movement or 
perspiration of a patient or the movement of an optical fiber that forms part of a source or 
detector. Thus, the results of an inverse calculation can include experimental systematic 
errors caused by measurement variations that are independent of the medium's properties 
of interest. The experimental systematic errors can also limit absolute spectroscopic 
measurements of the optical properties at a particular spatial location, i.e., absolute, rather 
than relative, values of absorption and scattering. 

[0009] Intemational Apphcation No. WO 01/19241 describes a calibration 

methodology for the diffuse optical measurements that corrects the transmittance 
measurements between a source and a detector for factors unrelated to sample properties. 
The calibration methodology is based on the same set of transmittance measurements that 
are subsequently corrected by the calibration, and then are used in imaging and/or 
spectroscopy applications. This calibration method involves a forward calculation for 
each of multiple source-detector pairs based on an approximate model of the sample, and 
a minimization of an expression that depends on the forward calculations and the 
transmittance measurements to determine self-consistent coupling coefficients for every 
source-detector pair. Once the coupling coefficients are determined, they are used to 
correct the transmittance measurements. An inverse calculation is performed on the 
corrected sample measurements to determine spatial variations in the optical properties of 
the sample. 
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SUMMARY OF THE INVENTION 

[0010] The present invention is preferably based on the fact that by allowing 

source and detector calibration factors to vary freely in the inverse calculation and be 
determined simultaneously with image reconstruction, image artifacts can be significantly 
reduced. 

[0011] In general, the present invention provides systems and methods for 

determining a distribution of one or more optical properties of a medium illuminated with 
radiation (e.g., electromagnetic radiation, infrared radiation, continuous-waive radiation 
etc.) from one or more sources, where the radiation exiting the medium is received by 
one or more detectors. The detectors may be adapted to obtain one or more parameters of 
the received radiation. Once the radiation is received and the parameters are obtained, 
one or more optical properties of the medium can be derived using at least one calibration 
factor and the received radiation, in which the calibration factors are variables. For a 
highly-scattering medium, the probability that photons entering the medium will scatter 
can significantly exceed the probability that photons entering the medium will be 
absorbed. Hence, employing diffusion approximation may be appropriate because the 
probability of scattering is at least an order of magnitude higher than the probability of 
absorption. Each source can be spatially separated from a respective detector, e.g. by two 
centimeters, and the separation selected can be limited by the experimental parameters 
(e.g., the number of sources and detectors, the radiation used, the size of the medium 
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etc.). Where the medium is human tissue, the separation may be one or more 
centimeters. 

[0012] The calibration factors used as variables in deriving one or more optical 

properties of the medium can include, for example, a source coupling factor, (which can 
be referred to interchangeably with a source strength or a detector gain), a detector 
coupling factor, a detector gain, a source location factor, and a detector location factor. 
[0013] The optical properties may be derived as a continuous-wave data and 

include, e.g., the absorption and scattering. However, the optical properties may also be 
derived, for example, as a frequency-domain or time-domain data. In case that the 
derived optical properties are in the frequency domain, each would have a respective 
amplitude and phase. Alternatively, if the derived optical properties are in the time 
domain, each would have a respective amplitude and temporal off-set. 
[0014] The or more of these optical properties can be derived by solving an 

inverse problem. The inverse problem can be solved by using, for example, a linear 
approximation model. One such exemplary linear approximation model used to derive 
the optical properties of the medium is a Bom Approximation. Another exemplary linear 
approximation model is a Rytov Approximation. In the Rytov Approximation model, the 
measured parameter (e.g., fluence), is defined as an exponential of a perturbation on the 
background fluence. A specific implementation of an embodiment of the method of the 
present invention using the Rytov approximation involves minimizing a least-squares 
expression for the difference between the logarithms of the theoretical and measured 
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fluence for each source-detector pair. The model can also scale the obtained parameters 
taken as input (e.g. fluence), to make those parameters dimensionless and to be on the 
same order as the calibration factors. 

[0015] Some examples of the property distributions obtained from performing the 

calculations on these measured parameters may include an absorption distribution, a 
scattering distribution, a distribution of both, etc. The determined property distribution 
can be displayed in one or more images, 

[0016] hi another embodiment of the present invention, a computer-readable 

medium containing a program is provided that can cause a processor to receive radiation 
scattered through an illuminated medium as input. Also, the processor may be 
programmed so as to obtain one or more parameters of the received radiation as input. 
The program can further cause the processor to derive one or more optical properties of 
the medium by using the input and one or more calibration factors. The program can 
subsequently cause the processor to determine a property distribution using the derived 
optical properties. The property distribution may be provided as output. The program 
can use a non-linear or Unear model to derive the optical properties of the media. Using 
the Rytov approximation model, the least-squares expression for the difference between 
the logarithms of the theoretical and measured fluence for each source-detector pair can 
be minimized so as to derive the optical properties, and determine the property 
distribution. The program can also scale all obtained parameters taken as input to make 
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them dimensionless and to be of the same order as the freely varying source and detector 
calibration factors. 



for determining the distribution of one or more optical properties of a medium 
illuminated with radiation from one or more sources. These systems have one or more 
radiation detectors for receiving radiation. These systems also include a processor for 
deriving one or more optical properties from the received radiation and one or more 
calibration factors. The processor also determines the distribution using one or more 
properties of the medium. 

[0018] The sources can have optical fibers and lasers, and the detectors may 

include optical fibers linked to photodetectors. To irradiate the medium, the system can 
use electromagnetic radiation, e.g., infrared radiation, or any continuous-wave radiation. 
The sources providing the radiation can be spatially separated from the detectors that 
receive the radiation. The spatially separated sources can extend in a plane opposite to a 
plane containing the detectors, and the spatial separation can be, e.g., two centimeters, 
and can be varied according to the experimental parameters (i.e., the number of sources 
and detectors, the radiation used, and the size of the medium). When the medium is 
human tissue, the separation can be one or more centimeters. The processor can also 
provide absorption coefficients, scattering coefficients, or both. 
[0019] In addition, the property distribution determined by the system of the 

present invention can be provided as an image on a display, e.g., a computer screen. The 
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In yet another embodiment of the present invention, systems are provided 
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model used to determine the property distribution can again be non-linear or linear. One 
such model includes the Rytov approximation and can involve, e.g., minimizing the 
square of the difference between the logarithms of the theoretical and measured fluence 
for each source-detector pair. The model can also utilize scaling parameters obtained 
from the signal as input to make them dimensionless and to be of the same order as the 
calibration factors. 

[0020] in still another embodiment of the present invention, systems are provided 

for determining a distribution of one or more optical properties of a medium illuminated 
with radiation from one or more sources. These systems have means for receiving 
radiation exiting the medium, and means for deriving one or more optical properties of 
the medium using one or more calibration factors and the received radiation, in which the 
calibration factors are variables. Means for determining a distribution in the medium 
using the derived optical properties are also provided. 

[0021] In yet another embodiment of the present invention, software systems are 

provided which, when executed on a processing device, may determine a distribution of 
one or more optical properties in a medium illuminated with radiation from one or more 
sources. These systems have a processing subsystem which, when executed on the 
processing device, configures the processing device to obtain parameters of received 
radiation exiting the medium, derive one or more optical properties of the medium using 
one or more calibration factors and obtained parameters, in which the calibration factors 
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are variables, and determine a distribution in the medium using one or more derived 
optical properties. 

[0022] Unless otherwise defined, all technical and scientific terms used herein 

have the same, or substantially similar, meaning as commonly understood by one of 
ordinary skill in the art to which the present invention belongs. Although methods and 
materials similar or equivalent to those described herein can be used in the practice or 
testing of the present invention, exemplary methods and materials are described below. 
In addition, the materials, methods, and examples are illustrative only and are in no way 
limiting. All cited references are incorporated herein by reference. 
[0023] In prior art systems and methods, several potential experimental 

systematic errors can hamper the ability to obtain images of the distribution of properties 
within a medium by applying radiation. While the characteristics of a system can be 
modeled prior to experimentation, slight fluctuations during experimentation, including 
for example changes in source strength or detector gain, can cause significant image 
artifacts. By including calibration in the image reconstruction algorithm, the present 
invention provides a substantial image improvement, despite the presence of possible 
considerable experimental systematic errors. One advantage of the methods and systems 
of the present invention is their ability to obtain diffuse optical tomography images 
having considerably reduced artifacts. The new imaging methods of the present 
invention can also be applied in electrical impedance tomography. These imaging 
methods can be utilized to image spatially varying optical properties of biological media. 
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e.g., tissue, with enhanced contrast. In addition, such novel methods can be applied to 
image non-biological media (e.g., plastics, ceramics, or liquids) to detect defects or 
impurities using radiation. Other features and advantages of the invention will be 
apparent from the following detailed description, and from the claims. According to the 
present invention, the source and detector calibration factors can be allowed to vary 
freely in the inverse calculation and be determined simultaneously with image 
reconstruction, and thus image artifacts can be significantly reduced. 

BRIEF DESCRIPTION OF THE DRAWINGS 
[0024] FIG. 1 is a top-level block diagram illustrating an exemplary embodiment 

of a method according to the present invention. 

[0025] FIG. 2 is a more detailed block diagram illustrating a step of deriving 

optical properties of the method shown in FIG. 1 . 

[0026] FIG. 3 is a schematic diagram of an exemplary embodiment of a diffuse 

optical tomography system according to the present invention. 

[0027] FIG. 4 is a schematic diagram of an exemplary source-detector pair for the 

system of FIG 1. 

[0028] FIG. 5 is a diagram of a first example of an experimental geometry 

derived using the present invention. 

[0029] FIGS. 6 A, 6B and 6C are exemplary images of simulated data obtained 

from the experimental geometry in FIG. 5 generated for 0%, 40% and 80% uncertainties, 
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respectively, using a normalized Rytov approximation and excluding the reconstruction 
of freely varying source and detector coupling factors. 

[0030] FIGS. 6D, 6E and 6F are exemplary images of simulated data obtained 

from the experimental geometry in FIG. 5 generated for 0%, 40% and 80% uncertainties, 
respectively, using a normalized Rytov approximation and including the reconstruction of 
freely varying source and detector coupling factors. 

[0031] FIGS. 7A, 7B, 7C and 7D are images of simulated data obtained from the 

experimental geometry in FIG. 5 generated for 0%, 40%, 80% and 120% uncertainties, 
respectively, using a Bom approximation and including the reconstruction of freely 
varying source and detector coupling factors. 

[0032] FIGS. 7E and 7F are images of simulated data obtained from the 

experimental geometry in FIG. 5 generated for 80% and 120% uncertainties, respectively, 
using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector coupling factors. 

[0033] FIGS. 8A, 8B, 8C, and 8D are images, similar to those of FIGS 7A-7F, of 

simulated data obtained from the experimental geometry in FIG. 5 generated for 10%, 
20%, 40% and 60% uncertainties, respectively, using a Rytov approximation with the 
reconstruction of freely varying source and detector coupling factors, but without the 
normalization. 

[0034] FIG. 9A is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 0.1 mm 
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using a normalized Rytov approximation and excluding the reconstruction of freely 
varying source and detector location factors. 

[0035] FIG. 9B is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 0.1 mm 
using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector location factors. 

[0036] FIG. 1 OA is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 0.5 mm 
using a normalized Rytov approximation and excluding the reconstruction of freely 
varying source and detector location factors. 

[0037] FIG. 1 OB is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for locaUzation uncertainty of e.g., 0.5 mm 
using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector location factors. 

[0038] FIG. 1 1 A is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 1.0 mm 
using a normalized Rytov approximation and excluding the reconstruction of freely 
varying source and detector location factors. 

[0039] FIG. 1 IB is a set of images of simulated data obtained from the 

experimental geometry in FIG. 5, generated for localization uncertainty of e.g., 1.0 mm 
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using a normalized Rytov approximation and including the reconstruction of freely 
varying source and detector location factors. 

[0040] FIG. 12 is an example of the experimental setup of the system according 

to the present invention. 

[0041] FIG. 13A is a set of images of data obtained from the experiment of FIGS. 

9A and 9B using a normalized Rytov approximation excluding the reconstruction of 
freely varying source and detector calibration factors. 

[0042] FIG. 13B is a set of images of data obtained from the experiment of FIGS. 

9A and 9B using a normalized Rytov approximation including the reconstruction of 
freely varying source and detector calibration factors. 

DETAILED DESCRIPTION 
[0043] The novel methods of the present invention calibrate sources and detectors 

as part of the inverse calculation for image reconstruction to reduce image artifacts. 

General Methodology 

[0044] Research in diffuse optical tomography has rapidly progressed in moving 

the imaging modality from the development stage of computer simulations and phantom 
experiments to the application stage of imaging animal and human subjects. However, 
little work has discussed the importance and solution of experimental details that degrade 
image quality including, without limitation, the treatment of boundary conditions for 
solving the forward problem, source and detector coupling to tissue, source and detector 
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positional uncertainties, and intrinsic tissue heterogeneity. The novel calibration 
techniques according to the present invention are based on the principle that, e.g., the 
inclusion of calibration factors as variables in the image reconstruction algorithm 
significantly improves the image quality. 

[0045] The calibration techniques of the present invention can be applied to 

determine a property distribution within a scattering medium, which is an ill-posed, non- 
linear inverse problem. A linear approximation can be obtained by assuming the 
variations may be described as small perturbations in absorption and scattering, 5]ia and 
b\Xs\ respectively, around known background values of absorption and scattering, \iao 
and Pso', respectively. As indicated previously, one or more parameters of diffuse light 
emitted through a scattering medium may be obtained and used as input in the inverse 
calculation. One of those parameters is preferably fluence. A linear approximation of the 
measured fluence 4> can be obtained as a perturbation ^^cn on a background fluence Oq. 
In imaging, the fluence <t> is measured experimentally and the perturbation cl^pen (which is 
caused by spatial variation in the radiative properties of the sample) should preferably be 
calculated. 

[0046] An accurate determination of the perturbation 4)pert prefers a calibration of 

a source and a detector, and a reasonable estimate of the background fluence <I>o. Errors 
in the determination of the caUbration factors are likely to result in a highly localized 
absorption and scattering perturbations appearing adjacent to the sources and detectors as 
a compensation for the errors in the image reconstruction process. This type of image 
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noise can be evidenced by high frequency spikes appearing preferentially near source and 
detector locations. 

[0047] The exemplary methods of the present invention allow for the calibration 

by enabling the calibration factors to vary freely in the inverse calculation, and 
determining them as part of the image reconstruction algorithm. There are various 
calibration factors, e.g., source and detector coupling factors s and d, respectively, and 
source and detector location factors Vs and rd , respectively. By considering the logarithm 
of the measured fluence O as input into the image reconstruction algorithm, the 
perturbation expert becomes linearly dependent on the logarithms of s and d. However, the 
perturbation Open can remain non-linearly dependent on the optical properties of the 
sample. To obtain a linear dependence on the optical properties, an approximation, e.g., 
the Rytov or Bom approximation, can be used, and the problem thus can become 
completely hnear. 

The Forward Problem: 

[0048] The radiative transport equation provides a rigorous theory to describe 

radiation emitted from a medium. It may be applied to the migration of photons through 
highly scattering media, such as biological tissue. This approach can indicate that near- 
infrared photons in highly scattering media e.g., human and animal tissue, essentially 
undergo a random walk, since the scattering probability greatly exceeds the absorption 
probability. The propagation of these photons through a highly scattering medium can 
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therefore be described by a diffusion approximation to the radiative transport equation. 
This exemplary diffusion can be provided as an equation as follows: 

-V.(D(r)VO(r,t)) + v^,(r)0(r,t)+^^J^ (1) 

dt 

where D(r) is the photon diffusion coefficient, 0(r,t) is the photon fluence at position r 
and time t, v is the speed of light in the medium, ]ia{r) is the absorption coefficient, and 
S(r,t) is the source distribution of photons (see, e.g., A. Ishimaru, Wave Propagation and 
Scattering in Random Media, Academic Press, Inc. San Diego, 1978; M. S. Patterson et 
al, Applied Optics 28:2331, 1989; R. C. Haskell et al, Journal of the Optical Society of 
America A, 11:2727, 1994). The photon diffusion coefficient Z)(r) is defined by the 
equation: 

D(r)=—--^ 7T' (2) 

3|ii/(r)+a iiijr) 

where laj'(r) is the reduced scattering coefficient. The coefficient a can be variously set 
to three, zero, or some other value (see, e.g., D. J. Durian, Optics Letters , 23:1502, 1998; 
R. Aronson and N. Comgold, Journal of the Optical Society of America A, 16:1066, 
1999). 

[0049] When focusing on the variations in the absorption, the value selected for a 

is moot, and therefore it can preferably be set to zero. 

[0050] When the distribution of properties within the medium varies spatially, 

two approaches can be used to find approximate solutions to Equation (1): the Bom 
approximation (see, e.g., A. Ishimaru, Wave Propagation and Scattering in Random 
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Media, Academic Press, Inc. San Diego, 1978; A. C. Kak and M. Slaney, Principles of 
Computerized Tomographic Imaging, IEEE Press, New York, 1988), and the Rytov 
approximation (see, e.g., A. C. Kak and M. Slaney, Principles of Computerized 
Tomographic Imaging, IEEE Press, New York, 1988). The Bom approximation can 
preferably be defined as: 

<^=^o+^^- (3) 
The Rytov approximation is preferably : 

<D=<D,exp(<Dp,rt). (4) 
[0051] Each approach decomposes the total fluence O into the background 

fluence Oo, which only depends on the background optical properties and \Xso\ and 
Opert, which is linearly related to the spatial variations in the optical properties 5|ia and 
Su^. For the Rytov approximation, assuming absorption variations only, the perturbation 
can be calculated as follows: 

<I>pen(r.>r,)=- / JcD^(r^,r)%^G(r,r,)rfr, (5) 

where and Yd are the respective positions of the source and detector, Do is the diffusion 
coefficient for the chosen a, and G is the Green's fiinction of the photon diffusion 
equation for the backgroimd optical properties provided by the boundary conditions. If 
the backgroxmd is homogeneous, G can be expressed analytically for some simple 
geometries (see, e.g., M. S. Patterson et al.. Applied Optics 28:2331, 1989; R C. Haskell 
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et al, Journal of the Optical Society of America A, 11:2727, 1994; S. R. Arridge et al, 
Physics in Medicine and Biology^ 37:1531, 1992). 

[0052] One exemplary geometry to which Equation (4) can be applied is a slab, 

i.e., a geometry where the sources are in a single plane on one side of the sample and the 
detectors are in a plane on the opposite side thereof. In this geometry, the extrapolated 
zero-boundary condition can be applied to the above-described radiative transport 
problem to solve for the Green's function G of the diffusion equation for the background 
optical properties (see, e.g., R. C. Haskell et al.. Journal of the Optical Society of 
America A, 11:2727, 1994). The background fluence Oo is related to the Green's function 
by the source and detector coupling factors. The source coupling factor s includes source 
power and a coupling to the medium. The detector coupling factor d includes detector 
gain and a coupUng to the medium. Using s and d respectively, to represent the source 
and detector coupling factors, the background fluence is as follows: 

^o(r.,r,)=5c/G(r„r,) (6) 

The Inverse Problem: 

[0053] By solving the inverse problem, the images of spatially varying optical 

properties can be generated from measurements of the fluence O. The inverse problem 
includes solving a least-squares equation. For the Rytov approximation, the expression to 
minimize can be as follows: 

FW=5[l"^7.e<,^,.- W - ln^A/...,f . (7) 
1=1 
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where the index / is summed over the measurements for each source-detector pair, <I>i(x) 

is provided by Equations (4) and (5) where x is a vector giving at each voxel 

position, and OMeas(x) is the measured fluence for each source-detector pair. 

[0054] For the case of fewer measurements than unknowns, the linear problem is 

underdetermined, and can be described by the (regularized) Moore-Penrose generalized 

inverse: 

x=-A''{AA^'^}jyy, (8) 
where lis the identity matrix, and X is the Tikhonov regularization parameter. Each 
element of the matrix A can be defined as: 



where r^., and r^j are the positions of the i"^ source and the i"' detector, respectively. The 
position of the j* voxel is expressed as Tj. Using the Rytov approximation, each element 
ofy can be as follows: 



y,=ln 



(10) 



[0055] In one exemplary approach, X can be set to 10*^ of the maximum eigen- 

value of /L4^(see, e.g., H. Dehghani et al., Physiological Measurement^ 20:87, 1999; V. 
Kolehmainen et al., Physiological Measurement, 18:289, 1997). 
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Image and Calibration Factor Reconstruction 

CALfflRATION FOR THE SOURCE AND DETECTOR 
COUPLING FACTORS AND IMAGE RECONSTRUCTION 

[0056] Image quality depends on accurate knowledge of the background optical 

properties, Uao and ]iso\ and the source and detector calibration factors s and d. 

Therefore, estimates of the matrices A and y similarly depend on these quantities. By 

considering \iao and ]Xso' as known properties, the inverse problem can be solved for the 

calibration factors. 

[0057] As described above, the measurements can include various errors caused 

by, for example, source and detector coupling to the medium, fluctuations in source 
power, variations in detector gain, and source and detector positional xmcertainties. 
Hence, the measurements can be calibrated for the source and detector coupling and 
intensity fluctuations (source and detector coupling factors), as well as for their positions 
( source and detector location factors). 

[0058] In the case when measurements should be caUbrated for the source and 

detector coupling factors s and d (the source and detector coupling and fluctuations in the 
source power and the detector gain), an augmented model of a particular measurement, 
may be written in terms of the unknown source coupling factors Sk, the detector coupling 
factors J/, and the absorption perturbations Suq/. The source and detector coupling 
factors can preferably refer to both coupling fluctuations and intensity fluctuations of the 
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sources and detectors. Using these variables, the particular measured parameter (in this 
case fluence) may be expressed as follows: 



3;,= In 



= hit,,,]+In[j,,,]+XAM.;. (11) 



where ^(0 and /(/) identify the source and detector, respectively, which are used for the i 
measurement. 

[0059] In case that numerous measurements are conducted, the measured 

parameter (fluence) can be written in a matrix form as: 

y=B^, whereB=[ASD], (12) 
where S is a source coupling factor matrix, D is a detector coupling factor matrix, and A 
is a matrix A rescaled by a factor of Uao. 

Factor ^ is defined by the formula: 



th 



In^; ••• Insj^^ Indj Inc/^^ 



(13) 



M'ao M'c 

where Nv is the number of voxels, Ns is the number of sources, and Nd is the number of 
detectors. 5]iaj can be scaled by to make the elements dimensionless and being of the 
same order as In 5 and In d. Scaling these quantities to make them dimensionless and of 
the same order as the logarithms of the calibration parameters yields a normalized 
formulation. 

[0060] As indicated previously, A is a rescaling of the matrix A whose terms are 

defined by Equation (9) by a factor of Uao. Thus, 
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(14) 



[0061] Matrices S and D have block diagonal form. For example, S has a one in 

the column for each measurement corresponding to source L and D has a one in 

column for each measurement corresponding to detector /. Thus, with four measurements 
between two sources and two detectors, the following matrix is defined: 



[SD]= 



10 10 
10 0 1 
0 110 
0 10 1 



CALroRATION FOR THE SOURCE AND DETECTOR LOCATION 
FACTORS AND IMAGE RECONSTRUCTION 

[0062] Conversely, if the measurements should be cahbrated for the source and 

detector locations, an augmented model of a particular measurement, yt, may be defined 

in terms of the uncertain source location 5rsk, uncertain detector location 5rdk and the 

absorption perturbations Siiaj- In this exemplary augmented model, yi may be expressed 

as: 



yi=\n. 



(15) 



where k(i) and /(/) identify the source and detector, respectively, used for the i"' 
measurement. It should be noted that 6rsk and 5rdi represent the differences from the 
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respective assumed values as used in and the real values for the measurement of In 
addition, 5rsk and 5rdi can be vector quantities representing the uncertainties in the x, y, 
and z coordinates of the corresponding source or detector. 

[0063] For this exemplary embodiment of the present invention, in case that 

numerous measurements are conducted, the measured parameter (fluence) can be written 
in a matrix form as follows: 

y=B^, where B=[APsPd] (16) 



where ^ is defined by: 



(17) 



in which Ny is the number of voxels, A^^ is the number of sources, and Nd is the number of 
detectors. Scaling 5]iaf by u^o and the positional unknowns by a makes the elements 
dimensionless. a can be selected so as to make the unknowns of the same order of 
magnitude, and in this manner, a is preferably an estimate of the expected error in the 
source and detector locations. Similarly, A is a rescaling of the matrix A whose terms are 
given by Equation (9) by a factor of iiao- Thus, 

A=^,,A. (18) 

[0064] Matrices Pg and Pd can have a block diagonal form. For example, Pg may 

be non-zero in the column for each measurement corresponding to source ky and Pd 
may be non-zero in the column for each measurement corresponding to detector /. 
These preferably have the same form as the [S D] matrix provided above with respect to 
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the image reconstruction including the calibration for the source and detector coupling 
factors. 

Frequency Domain Formulation 



using continuous-wave radiation. In the frequency domain, the measured parameters of 
amplitude and phase may also be obtained. If the measured parameters are in the 
frequency domain, then, for example, the source and detector coupling factors have 
corresponding unknown phases 6^ and 0^, respectively. They also likely have unknown 
amphtudes As and A^. The In s and In d, which appear in Equation (13), can be complex, 
as are the corresponding elements in the matrix A. To separate the real and imaginary 
portions of the complex matrix, the real portions can be stacked on top of the imaginary 
portions, i.e., 



[0065] 



The imaging method according to the present invention is not limited to 




(19) 



[0066] 



In the frequency domain, should be adjusted as follows: 



4= 



InA ••• ln^„ ••• G„ hiAj • • In^ 0. •■■ 9 



ao 



ao 



(20) 
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[0067] For the case of four measurements between two sources and two detectors, 

the complex matrix can be: 



[SD]= 



10/010/0 
10/010/0 
0 10/010/ 
0 10/010/ 



(21) 



[0068] The example of a formulation for the case where there are uncertain 

source and detector locations si substantially similar to the example provided above. In 
particular, Equation (19) would be provided as follows: 



B = 



re(A) re(P3) re(pj" 
im(A) im(P3) im(P^) 



(22) 



Time Domain Formulation 

[0069] The imaging method according to the present invention can also be 

applied to measurements taken in the time domain. To obtain time domain data, each 
source provides a temporally coherent pulse of radiation, e.g., a picosecond light pulse, 
and the detectors are time-gated to measure the temporal delay of the radiation in addition 
to its intensity. For a general reference on such time-domain techniques in diffuse optical 
tomography see, e.g., M.S. Patterson et a!.. Applied Optics 28:2331, 1989; S.R. Arridge, 
Inverse Problems, 15:841, 1999. Therefore, to apply the new calibration methods to 
time-domain data, the calibration factors s and d are described by corresponding 
unknown amplitudes and temporal off-sets. The unknown amplitudes in the time-domain 
are treated in the same fashion as in the continuous-wave domain. The unknown 
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temporal off-sets represent a non-linear perturbation and can be reconstructed using a 
non-linear form of the model. 

[0070] The example of a formulation for the case where there are uncertain 

source and detector locations si substantially similar to the example provided above. 

Exemplary Illustrations of the System and Method of the Present Invention 
[0071] Referring to FIG. 1, a top-level block diagram presenting an exemplary 

embodiment of a method according to the present invention is illustrated. Initially, the 
radiation which exits from the medium is received (step 1 10). Then, in step 120, one or 
more optical properties of the medium are derived using the radiation received in step 
110 and one or more calibration factors, wherein the calibration factors are variables. 
Subsequently, in step 130, a distribution of the optical properties in the medium is 
determined using the optical properties derived in step 120. The details of steps 110- 
130 are described in detail herein above. 

[0072] Further details for step 120 of deriving one or more optical properties of 

the medium are additionally illustrated in FIG. 2 and described below. In particular, 
parameters of radiation are obtained in step 210. One of such parameters can be fluence. 
The detectors with which the radiation exiting the medium is received may be adapted to 
directly measure the parameters of the radiation, e.g., fluence. Alternatively, the 
parameters may be calculated using the data corresponding to the received radiation. 
Once the parameters are obtained in step 210, the formulation of each obtained parameter 
can be augmented to include the calibration factors as variables (step 220). As described 
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above, the calibration factors may include the source coupling factor, detector coupling 
factor, source location factor and detector location factor. Subsequent to augmenting the 
formulations of the obtained parameters of radiation, one or more optical properties of the 
medixmi are derived as output by solving an inverse problem in which the augmented 
formulations are used as input (step 230). 

General Svstem 

[0073] FIG. 3 shows an exemplary embodiment of the system capable of 

executing the method according to the present invention in which experimental 
systematic errors may cause significant image artifacts. Such system can be a diffuse 
optical tomography ("DOT") system 300. The system 300 includes an array of spatially 
separated hght sources 310 and spatially separated detectors 320. During use, the array 
of sources 310 and detectors 320 can be positioned over a sample 350 to be imaged, e.g., 
a patient's head or breast. A controller 330 connected to the Hght sources 310 can 
sequentially trigger them to forward or provide light into the sample 350, so that a highly 
scattering medium may cause the light to become diffuse within the sample. For each 
sequentially triggered source, each detector 320 receives light that reaches it through 
sample 350 to obtain the measured fluence O. The controller 330 is also connected to the 
detectors 320, and selectively channels the signals from the detectors 320. An analyzer 
340 is connected to the controller 330 and analyzes the signals generated by detectors 
320. 
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[0074] For example, in one exemplary embodiment shown in FIG. 4, each source 

310 may include a diode laser 400 for producing optical radiation and a lens 420 for 
couphng the optical radiation into an optical fiber 410. The optical fiber 410 has an end 
415 adjacent to the sample 450 for directing the optical radiation into the sample 450. 
Each detector 420 includes an optical fiber 452 having an end 455 adjacent to the sample 
450 for receiving the optical radiation exiting the sample 450, a photodetector 460 for 
measuring the intensity of the optical radiation received by the optical fiber 452, and an 
amplifier 470 for amplifying the output of the photodetector 460. 
[0075] To obtain information about the optical properties of the sample 450 from 

the fluence, it is preferable to consider the measured fluence O as representing the 
perturbation $pert on a background fluence O©. According to Equation (4), this can be 
accomplished by using the Rytov approximation. By minimizing the least-squares 
expression as provided in Equation (7), the analyzer 340 can simuhaneously obtain the 
optical properties of the sample and conduct source and detector calibration. 
[0076] In other embodiments of the present invention, the light source can include 

a laser other than a diode laser, e.g., an ultrafast laser, or may be an incoherent soiu"ce. 
Also, the sources can include a common light source that selectively couples light into 
multiple fibers that deliver the light to spatially separated locations on the sample. 
Alternatively, the sources need not include optical fibers at all. For example, the lasers 
themselves can be positioned adjacent to the sample or can include beam delivery optics 
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to direct the light to the sample through free space. Furthermore, the light sources can 
provide light at multiple wavelengths by including, for example, multiple diode lasers. 
[0077] In the description provided below, it can be assumed that the sources 

provide continuous-wave optical radiation, and that the detectors measure the intensity of 
the optical radiation. 

[0078] It should be understood that the calibration techniques described herein 

can also be applied to other measurement techniques in which the sources do not provide 
continuous-wave radiation. For example, in some techniques, the amplitude of optical 
radiation provided by the sources can be modulated to create photon density waves in the 
sample, and the detectors are configured to measure the amplitude and phase of the 
photon density waves after propagation through the sample. For a general reference on 
"DOT" techniques with modulated optical radiation, see, e.g., M.A. O'Leary et al., Phys, 
Rev, Lett. 69, 2658 (1992). Furthermore, in other techniques, each source can provide a 
temporally coherent light pulse, for example, a picosecond pulse, and the detectors may 
be time-gated to measure the temporal delay of the difftise light pulse in addition to its 
intensity. For a general reference on such time-domain DOT techniques, see, e.g., M.S. 
Patterson et al., Appl. Opt, 28:2331, 1989, and S.R. Arridge, Inverse Problems, 15:841, 
1999. 

Exemplary Applications 

[0079] The methods and systems of the present invention provide the ability to 

obtain images of spatially varying optical properties within highly scattering media, such 
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as human and animal tissue, and can therefore be used to address many new biomedical 
problems and applications. For example, these novel methods and systems can be used to 
image a region of tissue to ascertain the presence of an early tumor, a small amount of 
bleeding, or an early aneurysm. In another exemplary application, multiple wavelengths 
can be used to evaluate local tissue concentrations of chemicals, such as hemoglobin, 
within tissue. Also, the methods and systems of the present invention can be further 
combined with use of a chemical stimulus to determine the chemical response. 
[0080] Non-biomedical appHcations for the methods and systems of the present 

invention also exist and can be used to image materials such as plastics, ceramics, and 
liquids for defects or impurities. In addition, the methods and systems of the present 
invention can be used for identifying regions of interest having distinct optical properties 
within highly scattering media, such as concealed objects, without invasive study of that 
media. For example, the methods of the present invention can be used to image an 
optical heterogeneity within a diffusive liquid. Furthermore, the novel calibration 
methods and systems can also be applied in electrical impedance tomography. 

Software Implementation 

[0081] The cahbration methods can be implemented in hardware, software, or a 

combination of both. The methods can be implemented via computer programs using 
standard programming techniques following the methods, equations, and Figures 
described herein. Program code can be executed using a processing arrangement which 
can receive input data to perform the functions described herein and generate output 
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information. The output information can be applied to one or more output devices, such 
as a display monitor, printer etc. 

[0082] Each program can preferably be implemented in a high-level procedural or 

object-oriented programming language to communicate with a general computer system. 
However, the programs can be implemented in assembly or machine language, if desired. 
In any case, the computer language can be a compiled or interpreted language. 
[0083] Each computer program can preferably be stored on a storage medium or 

device (e.g., ROM, magnetic diskette, or optical disc) readable by a general or special 
purpose programmable processing arrangement, for configuring and operating the 
computer when the storage media or device is read by the processing arrangement to 
perform the procedures and techniques described herein. The computer program can also 
reside in cache or main memory during program execution. The processing methods can 
also be implemented as a computer-readable storage mediimi, configured with a 
computer program, where the storage medium so configured can cause the processing 
arrangement to operate in a specific and predefined manner so as to perform the functions 
and techniques described herein. 

[0084] For example, referring again to FIG. 3, the analyzer 340 includes a 

processor 370, an input/output control card 360, a user interface 390 such as a keyboard 
and monitor, and memory 380. The memory 380 stores a program 385 specifying the 
steps of the calibration method. When executed, the program 385 causes the processor 
370 to carry out the steps of the method according to the present invention. 
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EXAMPLES 

[0085] The methods and systems of the present invention are further described in 

the following examples, which do not limit the scope of the invention described in the 
claims and are provided merely for illustrative purposes. 

Example 1 — Simulation of Image Reconstruction 

[0086] Several advantages of the methods of the present invention were 

demonstrated by simulations. The results obtained utilizing a normalized Rytov 
approximation without the reconstruction of calibration factors (i.e., without using the 
new calibration methods), a normalized Rytov approximation with reconstruction, a Bom 
approximation with reconstruction, and a Rytov approximation without normalization 
were compared. A diagram of a first example of an experimental geometry is shown in 
FIG. 5. 

[0087] The computer simulation involved the transmission through a 6-cm thick 

slab with background optical properties Uso' = 10 cm"^ and Uao, = 0.15 cm\ A 1.6-cm 
diameter absorbing object 550 with Uso' = 10 cm*' and iiao = 015 cm"' was centered at (x, 
y, z) = (1, -1, 3) cm. Sixteen sources 510 were programmed in a four-by-four grid at z = 
0, spaced 2 cm apart to span x = -3 to 3 cm and y = -3 to 3 cm. Similarly, sixteen 
detectors 520 were located in a four-by-four grid at z = 6, spanning x = -3 to 3 cm and y 
= -3 to 3 cm in 2-cm steps. All simulated measurements were made with continuous- 
wave sources. 
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[0088] The Full Bom expansion was used to solve the forward problem. A 1.6 x 

1.6 X 1.6 cm cube 555 divided into 7x7x7 voxels was centered over the 1.6-cm 
diameter absorbing object 550. 

[0089] Initially, the source and detector coupling factors were varied, namely, the 

source and detector amplitudes. The source and detector amplitudes were chosen 
randomly from a normal distribution with a mean of 1 . The effect of different standard 
deviations (0%, 40%, and 80%) was investigated. The normal distribution was biased by 
discarding negative amplitudes. A separate instance of the source and detector amplitude 
variation was considered for each standard deviation. There was no additive 
measurement noise (i.e., shot noise or electronic noise) in the simulated data, just the 
model error associated with the source and detector amplitudes. 
[0090] For the inverse problem, voxel centers were distributed from x = -3 to 3 

and y = -3 to 3 in 0.5-cm steps and z = 0.5 to 5.5 in 0.5-cm steps. There were 1859 
voxels, each with dimensions of 0.5 x 0.5 x 0.5 cm. The source and detector ampHtudes 
were assumed to be all equal to 1. Two different inverse problems were considered: (1) 
Rytov without coupling, and (2) Rytov with coupling. In each case, Tikhonov 
regularization was used to obtain the pseudo-inverse, since the system matrix is under- 
determined and ill-conditioned. 

[0091] FIGS. 6A - 6F show the results obtained with the normaUzed Rytov 

approximation both excluding the reconstruction of the source and detectors coupling 
factors (FIGS 6A-6C), and including such factors (FIGS. 6D-6F). FIGS. 6A and 6D 
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show results for 0% uncertainty in source and detector strengths, FIGS. 6B and 6E depict 
images for 40% uncertainty, and FIGS. 6C and 6F reflect 80% uncertainty. All images 
span X and Y from -3 to 3 cm, and Z-slices are vertically arranged from 0.5 (top) to 5.5 
cm (bottom). 

[0092] The results in FIGS. 6 A - 6C show images without the reconstruction of 

the source and detectors coupling factors, i.e., the effect of neglecting experimental 
systematic errors in the calibration. With no variance (i.e., no systematic error in the 
sources and detectors), the object is localized to the correct voxel in X and Y and is 
within 1 voxel (0.5 cm) in Z. The amplitude of the object is slightly underestimated due 
to blurring. While it is expected that 5iia = 0.10 cm\ the reconstruction yields 5iia = 
0.098 cm"\ The image contains artifacts near the sources and detectors. These artifacts 
disappear when the Bom approximation is used for the forward problem, indicating that it 
arises from mismatch between the exact forward solution and the Rytov approximation 
for the inverse model (i.e., it results from model noise unrelated to source and detector 
calibration factors). 

[0093] As the standard deviation in the calibration factors increases from FIG. 6A 

to FIG. 6C, it is still possible to ascertain the object, but this deviation has reduced the 
contrast relative to the surrounding voxels, and its amplitude is smaller than that in the 
voxels closest to the sources and detectors. The artifact amplitudes may greatly exceed 
the displayed contrast, and have been truncated to preserve contrast sensitivity sufficient 
to reveal the object of interest. The model noise caused by an incorrect calibration of the 
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source and detector calibration factors (i.e. systematic errors in s and d) clearly degrades 
the image quality by increasing the artifacts. 

[0094] FIGS. 6D - 6F show the reconstruction results where the source and 

detector coupUng factors are simultaneously determined with the optical properties of the 
sample using the novel calibration methods of the present invention. Without the 
variance in the calibration factors, the location of the object remains properly resolved in 
X and Y coordinates and is within 0.5 cm in Z coordinate. The resolution in FIGS. 6D - 
6F, with the simultaneous reconstruction of the calibration factors, is poorer than in 
FIGS. 6A - 6C, due to its lacking this simultaneous reconstruction. The presence of more 
unknowns in the inverse calculation using simultaneous reconstruction explains the 
decrease in the resolution. By contrast, the artifact shown in FIG. 6A due to the error 
between the exact and the Rytov approximation to the forward problem is reduced using 
the simultaneous reconstruction. The reconstruction appears to compensate for the model 
mismatch in the Rytov approximation. 

[0095] Comparing the images generated using the simultaneous reconstruction of 

the source and detector coupling factors (shown in FIGS. 6D - 6F) to the corresponding 
images generated without reconstruction (shown in FIGS. 6A - 6C) reveals a significantly 
improved image quality using the reconstruction. For both 40% and 80% variance, the 
image quality using the reconstruction in fact remains as good as the case where there is 
no variance. The minor differences are masked within the contrast sensitivity. While the 
reconstruction may not accurately determine the independent source and detector 
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coupling factors with the sample, it assists with the accurate determination of the source- 
detector coefficient products (e.g., Skdi) for each measurement. Despite the uncertainties 
of 80%, the results shown in FIG. 6D - 6F indicate that the source-detector coefficient 
products can be determined within an accuracy of a few percent. In addition, it has been 
observed that, when additive shot or electronic measurement noise is included in the 
simulation, there may be a greater sensitivity to the measurement noise as the source and 
detector variance increases. 

[0096] FIGS. 7 A - 7F show images of simulated data obtained from the 

experimental geometry in FIG. 5 generated for 0%, 40%, 80% and 120% uncertainties, 
respectively, using a Bom approximation including the reconstruction of freely varying 
source and detector coupling factors. FIGS. 7E and 7F are images of simulated data 
obtained from the experimental geometry in FIG. 5 generated for 80% and 120% 
uncertainties, respectively, using a normalized Rytov approximation including the 
reconstruction of freely varying source and detector coupling factors. FIGS. 7 A - 7F are 
useful to compare use of the Bom approximation with the reconstruction of the source 
and detector coupling factors to use of the normalized Rytov approximation with the 
reconstruction. FIGS. 7A - 7D show results for 0% uncertainty to 120% uncertainty in 
the source and detector strengths. By comparison, FIGS. 7E and 7F depict images for 
80% and 120% uncertainty using the Rytov approximation. The Z-slices are vertically 
arranged with each image spanning X and Y axes from -3 to 3 cm. 
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[0097] 



Images generated using the Bom approximation are clearly inferior to 



those generated using the Rytov approximation for the same level of uncertainty. Indeed, 
the images generated using the Rytov approximation for uncertainties of 80% (FIG. 7E) 
and 120% (FIG. 7F) resemble those generated using the Bom approximation at 0% (FIG. 
7A) and 40% uncertainty (FIG. 7B) more than images using the Bora approximation at a 
comparable level of uncertainty (FIGS. 7C and 7D). The standard Bom implementation 
of the image reconstmction problem is non-linear with respect to the source and detector 
coupling factors. Therefore, if the initial estimate for the coupling factors is off by more 
than approximately 10%, simultaneously solving for the coupling factors may not lead to 
an improved result. By using the Rytov implementation (i.e., taking the logarithm of the 
data), the dependence on the source and detector factors is made Unear, and it is then 
possible to accurately reconstruct the coupHng factors despite an initial guess, possibly 
off by several hundred percent. 

[0100] While merely substituting the Rytov approximation for the Bom 

approximation offers an improved image quality, the normalization of the different 
factors in the inverse calculation (i.e., absorption, scattering, and calibration factors) can 
also contribute to the enhanced image quality. The normalization has a beneficial effect 
because the reconstmction of several different factors simultaneously is difficult if the 
factors differ in magnitude by an appreciable amount. Using the normalization, the 
reconstmction can proceed accurately over a larger range of values. FIGS. 8A - 8D show 
the exemplary images using the Rytov approximation with the reconstmction of the 



NY02:36496M 



-38- 




A34927- 069225.01 10 



source and detector coupling factors, but without the normalization. FIG. 8A shows the 
results for 10% uncertainty, FIG. 8B for 20%, FIG. 8C for 40%, and FIG. 8D for 60%: 
FIGS. 8A - 8D have the same X-Y scale and arrangement of Z shoes as in FIGS. 7A to 
7F for an ease of comparison. Comparing FIGS. 8A - 8D with FIGS. 7E and 7F reveals 
that normalizing results in the improved image quality, even at higher levels of 
uncertainty. 

[0101] Alternatively, a different set of calibration factors may be varied. For 

example, the source and detector locations may be varied randomly with a normal 
distribution with a mean of 0. The effect of different distribution widths (0.1, 0.5, and 1.0 
mm) was investigated. A separate instance of the source and detector amplitude variation 
was considered for each distribution width. There was no additive measurement noise 
(i.e. a source or detector electronic noise) in the simulated data, and only the model error 
associated with the source and detector locations was added. 

[0102] For the inverse problem, voxel centers were distributed from x=-3 to 3 and 

y=-3 to 3 in 0.5 cm steps and z=0.5 to 5.5 in 0.5 cm steps. Approximately 1859 voxels 
were provided with dimensions of 0.5 x 0.5 x 0.5 cm. The vector ^ in Equation (14) was 
initialized to zeros. Two different inverse problems were considered: (1) the Rytov 
approximation without the locahzation, and (2) the Rytov approximation with the 
localization. In each case, the Tikhonov regularization was used to obtain the pseudo- 
inverse since the system matrix was under-determined and ill-conditioned. 
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[0103] FIGS. 9A-1 IB illustrate the results obtained with the Rytov approximation 

excluding and including the reconstruction of the coupling coefficients for optode 
localization errors of 0.1, 0.5 and 1.0 mm, respectively. Referring to FIGS. 9A and 9B, 
the images of the simulated data obtained fi-om the experimental geometry shown in FIG. 
5 using a normalized Rytov approximation excluding and including the reconstruction of 
freely varying source and detector location factors, are illustrated for localization 
uncertainty of 0.1 mm. With the localization error of 0.1 mm, it is possible to locate the 
position of the absorbing heterogeneity despite the presence of image speckle near the top 
and bottom surfaces of the scattering medium. 

[0104] With the localization errors of 0.5 mm and greater (See FIGS. lOA and 

lOB, and FIGS. 1 1 A and 1 IB, the surface speckle noise obscures the presence of the 
absorbing heterogeneity. The speckle noise may arise directly firom the localization error. 
In most, if not all, cases, including the calibration of optode location in the inverse 
problem significantly reduces the surface speckle noise revealing the absorbing object. 

Example 2 — Phantom Results 

[0105] The actual image results of a phantom have also been obtained, and FIG. 

12 shows the experimental setup for imaging the phantom. A phantom box 1210 has a 
set of sources and a set of detectors on its top and bottom plates, respectively. The 
phantom box 1210 is connected with an RF Imager 1230 with a set of detector fibers 
1220. The phantom box 1210 is also connected with a controller 1280 via a set of source 
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fibers 1270, A processing arrangement 1290 is coupled with the controller 1280. The 
phantom box 1210 is further coupled with a beaker 1240 and a pump 1250. 
[0106] The phantom box 1210 was filled with 0.5% intralipid solution, and a 2- 

cm diameter heterogeneity sphere was embedded in this solution. The content of the 
sphere was controlled by the circulation channel outside the phantom box 1210, which 
included the pump 1250 and the beaker 1240 for ink titration. In the experiment, a slab 
geometry was formed, with the slab having a thickness of 5.1 1 cm. The center of the ball 
was located 2.35 cm firom the bottom of the slab. There were sources on the top plate (z 
= 0 cm) and detectors on the bottom plate (z = 5.1 1 cm). The signal to the sources was 
deUvered via the source fibers 1270. The signal was Hght with a wavelength of 830 nm 
and was modulated at 70 MHz. The amplitude and phase of the diffusely transmitted 
light were measured. Data was collected firom the detectors via the detector fibers 1220 
into the radio firequency imager 1230. The processing arrangement 1290 processed the 
results and generated the images. 

[0107] The images were first generated using the Rytov approximation both 

excluding and including the reconstruction of the source and detector coupling factors. 
Images were median filtered using standard techniques because this was found effective 
in removing speckle artifacts. FIG. 13A shows the images obtained excluding the 
reconstruction of the source and detector coupling factors, and FIG. 13B shows the 
images generated including the reconstruction. The X and Y ranges are in centimeters, 
and Z-slices are vertically arranged fi-om 0.01 (top left) to 5.501 cm (bottom right) in 
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0.50-cm steps. It is observed that including reconstruction of the source and detector 
coupUng factors improves image quality. In particular, the detail of the heterogeneity 
sphere can be seen in FIG. 13B for Z=2.01 to Z=3.51 cm, but it may not be observable 
over the corresponding Z-range in FIG. ISA. 

[0108] It is to be understood that while the invention has been described in 

conjunction with the detailed description thereof, the foregoing description is intended to 
illustrate and not limit the scope of the invention, which is defined by the scope of the 
appended claims. Other aspects, advantages, and modifications are within the scope of 
the following claims. 
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